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Abstract — An efficient one dimensional 3-noded finite element 
model has been developed for the vibrational analysis of 
composite beam for various boundary conditions, using the 
efficient layerwise zigzag theory. The results are compared with 
2-noded beam element model. To meet the convergence 
requirements for the weak integral formulation, fifth power 
Hermite interpolation is used for the transverse displacement and 
quadratic interpolation is used for the axial displacement and 
shear rotation. Each node of an element has four degrees of 
freedom. The formulation is validated by comparing the results 
of the 2D finite element (2D-FE) for the simply supported beam. 
The present zigzag finite element results for natural frequencies 
and mode shapes of the beam are obtained with one-dimensional 
finite element (1DFE) codes developed in MATLAB. This 
comparison establishes the accuracy of zigzag finite element 
analysis for natural frequencies of the symmetric laminated 
beams, and helps to obtain natural frequencies of fixed beam and 
cantilever beam with less error. 

Index Terms — composite laminates, sandwich beam, 
Vibration analysis, zigzag theory, FEM, MATEAB. 


I. INTRODUCTION 


Composite structures are increasingly used in areas like 
automotive engineering and other applications as they posses 
lower weight and higher strength and stiffness than those 
composed of other metallic materials. For design of composite 
and sandwich beams accurate knowledge of deflection and 
stress under static and dynamic loadings, natural frequencies, 
mode shapes are required. Exact elasticity solutions [1-3] have 
been provided for static, free and forced vibration cases. 
Discrete layer theories with layer wise displacement 
approximation are quite accurate but computationally 
expensive as the number of basic variables depends on number 
of layers. 
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Among most discrete layer theories [4, 5] the shear stress 
continuity at layer interfaces, is violated. Kapuria et. al. [6, 11] 
presented an assessment of zigzag theory for laminated 
composite beams by giving analytical solution for simply 
supported end conditions only Benjeddou [7] has presented 
finite element modeling of adaptive structures. Kapuria et. al. 
[8] presented a novel finite element model of efficient zigzag 
theory for static analysis of hybrid piezoelectric beams. They 
presented finite element analysis of hybrid piezoelectric beams 
under static electromechanical load using zigzag theory. They 
also compared then results with 2D finite element results 
obtained using ABAQUS to establish the accuracy zigzag 
theory. Navier type solutions for simply supported beams were 
presented in refs. [6, 9] which did not provide finite element 
formulation of zigzag theory. One of the authors [10] recently 
presented efficient layer wise finite element model for 
dynamic analysis of laminated piezoelectric beams. Two 
noded beam element model developed by Alam et. al [13] is 
an initiative approach in the field of F E analysis. 

This work considers a three nodded finite element model 
for dynamic analysis for composite beams based on zigzag 
theory 6 in which the shear traction condition at the top and 
bottom and the transverse shear continuity condition at the 
layer interfaces are satisfied. 5 th Hermite interpolation function 
[12] is used for deflection and quadratic interpolation is used 
for the axial displacement and rotation. The finite element 
formulation and the MATLAB code developed thereof are 
validated by comparison of the results with the results 
obtained using ABAQUS software. The one-dimensional 
finite element (1D-FE) results for propped beam are compared 
with 2D finite element (2D-FE) results. 


II. FINITE ELEMENT MODEL USING 3-NODED 
BEAM ELEMENT 

With reference to two noded element model [13], three nodded 
elements are used for the displacement variables. The primary 
variables, u 0 wo, Wo, within an element are expressed in terms 
of their nodal values using appropriate polynomial 
interpolation functions. The highest derivatives of u 0 , Wo, Wo„ 
appearing in the Variational equation [13] are u 0x , Wo, X x, Wo.x- 
To meet the convergence requirements of the finite element 
method, u 0 , wo lX , Wo, must be continuous at the element 
boundaries. Hence w 0 is expanded using 5 th power 
interpolation along x in terms of the nodal values of w 0 , and 
w 0 ,x ■ Similarly, quadratic Lagrange interpolation along x is 
used for u 0 and Wo in terms of their nodal values. Thus at the 
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element level, each node will have four degrees of freedom u 0 , 
Wo, Wo, x , Wo, for the displacements. The values of an entity (...) 
at the nodes 1, 2 and 3 are donated by (...)i, (...) 2 , and (...) 3 
respectively. 


Le 

T e = \sU eT [Bl.IB m U e +B T D.BU e -B T m2 F 2 ]&c =0 

o 

(13) 


The following interpolations of u 0 , Wo, Wo, have been used in 
terms of the nodal values and the shape function matrices N 

and N : 
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N )=2 x 2 /hf -x/L e , 
Ni= x/L P +2 x 2 /L p 2 


N 2 =l- 4(x/L e ) 2 , 


Similarly, other shape functions are derived using 
interpolation [12] as: 

N ,= x 2 /L 2 -5/4.x 3 / L e 3 -l/2.x 4 /L e 4 +3/4.x 5 / L e 5 

N 2 = 1/4 ,x 2 /L e +1/4. x 3 /L e 2 -l/4.x 4 /L 3 +l/4.x 5 / L e 
( 12 ) 


N„ M x , P„ and Q x are substituted from ref. [13] to obtain 
general equation after integration as: 

M e U e + K e U e = P e (14) 

where U e =[U e Q T Wq T l// e Q ] T , M e = J B T jB m dx . , 

-Le 

P e = [ 0 N t F 2 0] , K' = ^B j IB dx. 

-Le 

The natural frequency for free vibration can be obtained from 
the above equation by making it aneigen value problem as 

Kflf- CO 2 M e U e =0 (15) 


The beam generalized displacement are related to 
displacement vector 

U (displacement vector) = [U Q w Q y/ Q ] 


+ =B„,U- ”‘ l u^B. 2 U 



0 



~N 

0 

0" 

~N 

0 " 






0 


B ,nl 


0 

-N x 

0 


-N x 

and B m = 


— 


0 


N 

m 



0 

0 

N 



\_ ml J 







0 

N 

0 


Refer appendix for the equations and matrices. 


III. RESULT AND DISCUSSION 


N 3 = I -2x/ L e 2 + x 4 /L, 4 , 

N 4 =x -2x 3 / L e / L e 4 
N 5 = x 2 / L e 2 +5/4. x 3 / L e 3 -l/2.x 4 /L e 4 -3/5.x 5 / L e 5 
N 6 = -l/4.x 2 / L e -1/4. x 3 / L e 2 +l/4.x 4 /L e 3 +l/4.x 5 / L e 4 

The integrand in the variational equation can be expressed 
as, [ref. 13, equaion no. 23] 


A highly inhomogeneous symmetric beam is analyzed for 
simply supported boundary condition. The stacking order is 
mentioned from the bottom. The beams is composite beam of 
material [11] with Y/ = 181GPa and Y 2 = Y, =10.3 GPa and 
Y 12 = V i3=-25 , V 2 j=.33 .consisting of four plies of equal 
thickness 0.25h. It has symmetric lay-up [0/90/90/0]. The 
density of materials of the beam is 1578 kg/m 3 [1 1]. Flexural 
natural frequencies and mode shapes of laminated composite 
beams have been computed by developing 1D-FE MATLAB 
program [11]. [Table 1] 
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Fig. 2 % error in CO , for 1 st longitudinal mode using 2 noded Fig- 3 % error in CO i for 1 st longitudinal mode using 3-noded 
beam elements beam elements 


Table 1 % error of natural frequencies of symmetric beam (b) 
simply supported boundary condition using FEM and 2D exact 
values 


n 

m 

s 

2D 

exact 

(Availbl) 

% error 2- 
noded 

% error 
3-noded 

i 

i 

5 

6.8060 

0.13 

0.1 


i 

10 

9.3434 

0.068 

0.04 


i 

20 

10.640 

0.056 

0.03 


i 

100 

11.193 

0.032 

0.02 

2 

i 

5 

16.515 

1.33 

.91 


i 

10 

27.224 

.14 

0.08 


i 

20 

37.374 

.047 

0.02 


i 

100 

44.477 

0.02 

.011 

3 

i 

5 

26.688 

4.41 

2.81 


i 

10 

46.419 

0.54 

0.31 


i 

20 

71.744 

0.022 

0.02 


i 

100 

98.988 

0.001 

0.001 


IV. CONCLUSION 


The present FE model is developed for vibration analysis of 
symmetric beam for various end conditions. The presented 
model can be used for computing frequencies for various end 
conditions. The comparison shows that the 1DFE model of 
zigzag theory yields very accurate results for mode shapes for 
simply supported beams. 
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V. APPENDIX 
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the constants are given by 
L e 

C 8 = \ N J- N J X ’ C 9 
—le 


L e 

\Nj.N„dx. 

—Le 


L e 

C 10 = jAxv •N'X'dx 

-Le 


c 8 = [ 38/ (3*L e ) , -64/ (3*L e ) , 26/(3*L e )] 

[ -64/ (3*L e ) , 128/ (3*L e ) , -64/(3*L e )] 
[ 26/ (3*L e ) , -64/ (3*L e ) , 38/(3*L e )] 


c 9 = [ 4/L e A 2, 13/L e , 0, 0, -4/L e A 2, 3/L e ] 

[-8/L e A 2, -24/L e , 0, 0, 8/L e A 2, -8/L e ] 
[ 4/L e A 2, ll/L e , 0, 0, -4/L e A 2, 5/L e ] 


[ 1273/ (70*Le A 3) , 
[ 779/ (70*Le A 2) , 
[ -64/ (5*Le A 3) , 
[ 96/ (7*Le A 2) , 
[ -377/ (70*Le A 3) , 
[ 121/ (70*Le A 2) , 


C10 — 

77 9/ (70*Le A 2) , 
58 6/ ( 35*Le ) , 
-32/ (5*Le A 2) , 
32/ (7*Le) , 
-331/ ( 7 0*Le A 2 ) , 
124/ (35*Le) , 


-64/ (5*Le A 3) , 
-32/ (5*Le A 2) , 
128/ (5*Le A 3) , 

0 , 

-64/ (5*Le A 3) , 
32/ (5*Le A 2) , 


96/ (7*Le A 2) , 
32/ (7*Le) , 

0, 

128/ (7*Le) , 
-96/ (7*Le A 2) , 
32/ (7*Le) , 


-377/ (70*Le A 3) , 121/ (70*Le A 2) ] 
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-331/ (70*Le A 2) , 
-64/ (5*Le A 3) , 
-96/ (7*Le A 2) , 
1273/ (70*Le A 3) , 
-569/ (70*Le A 2) , 


124/ (35*Le) ] 
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[ (14*L e )/15, 
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- (67*Le A 3) /6930, 
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